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Abstract 

Nonlinear two-point boundary value problems arise in numerous ar- 
eas of application. The existence and number of solutions for various 
cases has been studied from a theoretical standpoint. These results 
generally rely upon growth conditions of the nonlinearity. However, 
in general, one cannot forecast how many solutions a boundary value 
problem may possess or even determine the existence of a solution. 
In recent years numerical continuation methods have been developed 
which permit the numerical approximation of all complex solutions of 
systems of polynomial equations. In this paper, numerical continu- 
ation methods are adapted to numerically calculate the solutions of 
finite difference discretizations of nonlinear two-point boundary value 
problems. The approach taken here is to perform a homotopy defor- 
mation to successively refine discretizations. In this way additional 
new solutions on finer meshes are obtained from solutions on coarser 
meshes. The complicating issue which the complex polynomial sys- 
tem setting introduces is that the number of solutions grows with the 
number of mesh points of the discretization. To counter this, the use 
of filters to limit the number of paths to be followed at each stage is 
considered. 
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1 Introduction 



Consider a two-point boundary value problem on the interval [a, b] C R, 

y" = f(x,y,y r ), (1) 

with boundary conditions y{a) = a and y{b) = fi. The standard central 
difference approximation with a uniform mesh may be used to approximate 
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solutions to In particular, let N be a positive integer, h := ^^j, and 
Xi := a + ih for i = 0, . . . , N + 1. Setting yo = « and yN+i = /?, the 
discretization of Q takes the form of the system T>^\ 

yo - 2yi + y2 =h 2 f(x uyu ^) 

y N -i - 2y N + y N+1 = h 2 f(x N , y N , m± ^^ ) 

A solution y{x) of (JJ) may then be approximated by an iV-tuple of real 
numbers (yi, . . . , yjv) such that yt « y(xj) for i = 1, . . . , N, 

Depending upon the nonlinearity /, equation may have no solu- 
tions, a unique solution, multiple solutions, or even infinitely many solu- 
tions. There are many existence theorems for solutions of such equations 
subject to growth conditions on /, but even when existence is known, the 
number of solutions often is not. Furthermore, a discretization such as T>n 
may have spurious solutions that do not converge to a solution to Q as 
N — > oo. On the other hand, if / is sufficiently smooth, a solution y to 

is eventually approximated with 0(h 2 ) accuracy on the mesh by some 
solution y € Mr. 

The purpose of the present paper is to give a relatively secure numerical 
technique for finding the solutions of a general class of two-point boundary 
value problems without requiring highly refined meshes. The technique in- 
volves performing successive homotopy deformations between discretizations 
with increasingly many mesh points, as suggested in 2 . By restricting our 
attention to problems having polynomial nonlinearity, including the case of 
a polynomial approximation to a smooth nonlinearity, we can often assure 
that all solutions are found at each stage of the algorithm. Even when we do 
not guarantee all solutions, our method generates multiple solutions that in 
test cases include approximations to all known solutions. While Grobner ba- 
sis methods (see jl]) or cellular exclusion methods (see [H]) could be applied 
to solve the polynomial discretizations, we chose to use homotopy contin- 
uation due to its ability to handle polynomial systems in many variables 
and the ease with which it allows us to generate solutions on a refined mesh 
from the solutions on the previous mesh. Although other numerical meth- 
ods treating two-point boundary value problems have been developed (see 
[7] and jlUj). such methods require satisfactory initial solution estimates. 
The present technique provides such initial estimates. 

Here is a sketch of our bootstrapping process, which will be discussed in 
more detail in the subsequent section: 

1. Find all solutions of the discretization T>n for some small N. The size 
of N needs only to be large enough that the discretization is consistent; 
it could be as small as N = 1. 
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2. Discard all unreasonable solutions, e.g., solutions which do not possess 
properties which exact solutions may be known to have. Let us denote 
the set of solutions which are kept by Vjv. 

3. If the mesh size is not yet sufficiently small or the cardinality of Vn 
has not yet stabilized, add a mesh point to obtain the discretization 
T>n+i- Use the solutions in Vn to generate solutions Vn+i of T>n+i 
and then return to Step 2. 

4. Once the mesh size h is sufficiently small and the cardinality of Vat 
becomes stable, refine the solutions with a fast nonlinear solver, using 
starting values obtained by interpolating the solutions in Vn- 

This paper focuses primarily on the implementation of Step 3 of the above 
scheme. In particular, we consider the homotopy function 



H N +i(yi, • • .,yN+i,t) := 

yo - 2yi + y2 

y N ^ 2 - 2yAT_i + yN 
y N -i - 2y N + Y N+1 (t) 
y N - 2y N+1 + Y N+2 (t) 



h(t) 2 f(x 1 (t),y 1 , 



y-2-yu 
2h(t) 



- h(t) 2 f[x N -i(t),y N . 



1; 



VN-yN-2 
2h(t) 



h{tff[x N {t),y N , YN+ ^ N - 1 
h(t) 2 f{x N+1 (t),y N+1 / N ^ yN ) 



with 



yo 

h(t) 

Y N+ i(t) 
Y N+ 2(t) 



= a 



(I - t)y N+ i + Pt 
0(1 -t) 



Xi(t) := a + ih(t), i = 1, . . . , iV + 1. 

At t = this is the system T>n+i- At t = 1, it can be interpreted as the 
system T>n with a new mesh point having the value yN+i at x = b and 
a new right-hand boundary at x = b + h(l) having value Y/v +2 (l)- The 
incompatibility of the old boundary condition at x = b and the new one 
at x = b + h(l) is accommodated by the presence of both yN+i and Y/v+i, 
which are not necessarily equal. As t goes from 1 to 0, the mesh points 
are squeezed back inside the interval [a, b], and the right-hand boundary 
condition y(b) = (3 is transferred from Y/v+i to Y/v+2 as Yn+i is forced to 
equal yjv+i. 

To find solutions of T>n+i, we use continuation to track solutions of Hn+i 
as t goes from 1 to 0. At t = 1, we have a list Vn of solutions (yi, . . . , yjv) 
satisfying the first N equations of Hn+i, while the final equation is 



y N - 2y N +i = h(t) 2 f b,y N+1 , 



-yN 
2h 
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which is the only place where un+i appears. For each solution of (y\, . . . , yj\[) 
in Vat, we may use this equation to find corresponding solution values for 
yN+i- These are the start points of continuation paths leading to solutions 
of V N+1 . 

The framework above does not change in any of its essentials if we pre- 
scribe a different function for Y/v+2(^)- For example, the constant function 
YN+2(t) '■= (3 was used for all examples below. Although other alternatives 
are theoretically feasible, none were tested. The essential feature of Yjv+2 (t) 
is that it goes to j3 as t goes to 0. 

By the implicit function theorem, a nonsingular solution y = y* to 
fljv+i(y 5 1) = will continue uniquely in the neighborhood of t = 1 to a 
nonsingular solution path y(t) satisfying H^ + \{y{t),t) = with y(X) = y*. 
This does not mean, however, that the path remains nonsingular all the way 
to t = 0, which is what we require to follow the path reliably with numerical 
continuation. To skirt this difficulty, as discussed in Chapter 7 of it is 
sufficient to insert a random 7 € C into the homotopy to obtain the variant 

Hn+i(ui, ■ ■ -,yN+i,t) := 

T(t)(yo-2y 1+ y 2 ) - h{tf f (*i(t), y u ^f) 

T(t) (yjv-2 - 2?Mr_i + y N ) 
r(t) {y N . 1 -2y N )+Y N+1 {t) 
r(t) {y N ~ 2y N+1 + 0) 

with 

r(t) := 7 2 t + (l-t) 

hit) := 7t (fe) + (1 " *) (M) 
Fiv + i(i) := (1 -t)yjv+i +^ 2 f3t 

Xi (t) := a + ih(t), i = l,...,N + l. 

The work of the second author was supported by the National Science 
Foundation under Grant No. 0105653 and Grant No. 0410047; and a fel- 
lowship from the Arthur J. Schmitt Foundation. The work of the third 
author was supported by the National Science Foundation under Grant No. 
0105653 and Grant No. 0410047; and the Duncan Chair of the University of 
Notre Dame. The work of the fourth author was supported by the National 
Science Foundation under Grant No. 0410047. 



- h (tff{x N ^{t),y N ^ y -^^) 
h(t) 2 f(x N+l (t),y N+1 ,^') 

(2) 



4 



2 The case of polynomial nonlinearity 

Let's specialize the homotopy of equation to the case when f(x,y,y') 
is a real polynomial p(y). Then, the right-most term of the i th entry in 
iJjv+i(y, t) becomes just h 2 (t)p(yi). This restriction to the polynomial case 
allows us to conveniently obtain the start points for HN+i(yi, . . . , yN+i, 1) = 
by solving the polynomial 



for yN+i given y^ from the solutions in Vn- 

Let d = degp(y). We see that, in general, over the complex numbers, 
we will obtain d values of yN+i for every point in Vn- Suppose that at each 
stage of the algorithm these all continue to finite, nonsingular solutions of 
T>n+i- Then, the solution list Vn will have d N entries. While this gives 
an exhaustive enumeration of the solutions of the discretized problem, the 
exponential growth in the length of the solution list cannot be practically 
sustained as N increases. However, it is often the case that most of the 
solutions at a given stage do not exhibit various properties required of solu- 
tions to the two-point boundary value problem at hand, leading to filtering 
rules. Depending upon the problem at hand, there are a variety of filtering 
rules that may be implemented to determine which solutions in Vn may be 
discarded as start solutions for the subsequent homotopy. 

For small N we can contemplate retaining all solutions. It is reasonable 
to ask whether the above procedure is guaranteed to generate all solutions of 
the discretized system. The answer, in general, is no, but we can say that if 
Vn has d N distinct, nonsingular solutions, then it is clear that all solutions 
have been found, as Bezout's theorem states that this is the greatest number 
possible. Indeed, for our test problems, we have found that this behavior is 
typical. 

For larger N, a filter becomes necessary. One that is always available 
is to take the discretization of the derivative of y'" = p'(y)y' of y" = p(y), 
and throw away y € Vn for which this is large. To get the discretization we 
could use the central difference approximations 



applied only at the mesh points y2, ■ ■ ■ , 2/jv-l- So we would throw away the 
point z = (yi, ... ,y N ) G Vn if 





Ui+i ~ Vi-l 
2h 



and 




Ui+2 ~ tyi+i + %-i ~ Vi-2 
2h 3 



N-l 




1=2 
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for some €2 > 0. Naturally, one drawback to such a filter is the need to 
specify €2- 

Other filters may be derived from known properties of the solutions of 
the problem at hand. For example, it may be known that solutions are 
symmetric about always positive, oscillate with a specific 

period, or exhibit some other easily-detected behavior. For example, a filter 
based on symmetry is considered in Section 3.2. Although one may be 
tempted to discard solutions having nonzero complex part, this is not a 
valid filtering rule. The problem in Section 3.3 below has non-real solutions 
in Vn that are tracked to real solutions in Vn+i- Similarly, it is possible that 
oscillating solutions may arise from a sequence of non-oscillating solutions 
and that similar problems may occur with other filters. Thus, the use of 
filters may be computationally beneficial, but with it comes the risk of not 
finding all real solutions to the problem. 

Thus we have the final version of the algorithm: 
Algorithm 1 

1. For N = 1, H± (yi,l) is a single polynomial in yi, which 
may be solved with any one-variable method to produce 
Vi. 

2. For N = 2, 3, until some desired behavior has occurred: 

(a) Form the homotopy Hn(vi, ■ ■ ■ , Vn, t). 

(b) Solve the last polynomial of Hn(vi, . . . , yN, t) for yN 
using each solution in Vjv-i, thereby forming the set S 
of the start solutions for Hn(vi, ■ ■ ■ ,VN,t)- 

(c) Track all paths beginning at points in S at t = 1. The 
set of endpoints of these paths is Vn- 

(d) If desired, apply a filter to Vn to reduce the number of 
paths to be tracked in stage N + 1. 

3. Refine the solutions with a nonlinear solver, if desired. 

3 Numerical experiments 

The following experiments were run using Bertini, a software package under 
development by the last three authors for the study of numerical algebraic 
geometry. Although Bertini was written to make use of multiprecision adap- 
tively, each of the following experiments ran successfully using only 16 digits 
of precision. 

In the following, N denotes the number of mesh points, SOLS(iV) de- 
notes the total number of solutions (real or complex), and REAL(iV) denotes 
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the number of real solutions. For N > 1, the number of paths tracked from 
stage N — 1 is d ■ SOLS(iV — 1). A solution is considered to be real if the 
imaginary part at each mesh point is zero to at least eight digits. 

3.1 A basic example 

As a first example, consider the following two-point boundary value problem 

y" = 2y 3 (3) 

with boundary conditions y(0) = \ arm = \- 

There is a unique solution, y = to ©• Our method produces 

one real solution among a total of 3^ solutions found for TV" = 1,...,9. 
Furthermore, the error between the computed solution and the unique exact 
solution behaves as 0{h 2 ). Refer to Table 1 for details. 



N 


Maximal error at any mesh point 


K 1 


Maximal error/ h' 2 


3 


1.570846e-04 


4.000000e-02 


3.927115e-03 


4 


1.042635e-04 


2.777778e-02 


3.753486e-03 


5 


7.069710e-05 


2.040816e-02 


3.464158e-03 


6 


5.348790e-05 


1.562500e-02 


3.423226e-03 


7 


4.078910e-05 


1.234568e-02 


3.303917e-03 


8 


3.230130e-05 


1.000000e-02 


3.230130e-03 


9 


2.624560e-05 


8.264463e-03 


3.175718e-03 



Table 1: Evidence of 0(h 2 ) convergence for Problem 



3.2 A more sophisticated example 

Consider the problem 

y» = -\(l + y 2 ) (4) 

with zero boundary conditions, y(0) = and y{l) = 0, and A > 0. 

According to , any solutions to this problem must be symmetric about 
x = ^, so we have a special filter. Furthermore, it is known that there are 
two solutions if A < 4, a unique solution if A = 4, and no solutions if A > 4. 
Without using a filter, the expected number of real solutions in the first 
and last cases were confirmed computationally (for A = 2 and A = 6 with 
N < 17), and the computed solutions were symmetric as anticipated. From 
Bezout's theorem, one would expect to find at most 2^ complex solutions 
at each stage N, and this is precisely the total number of complex solutions 
found. When A ~ 4, the Jacobian of the associated polynomial system is 
rank-deficient, so regular path-tracking techniques fail. 



7 



0.2 0.4 0.6 0.8 1 

x 



Figure 1: The real solutions of © with N = 20. 

Tracking all 2 17 paths for N = 17 took just under an hour of CPU time 
on a single processor Pentium 4, 3 GHz machine running Linux. At this 
rate, ignoring the time-consuming data management part of the algorithm, 
it would take well over one year to track all 2 30 paths for N = 30 mesh 
points. As discussed in Section 2, filtering rules may be used to dramatically 
reduce the number of paths to be tracked at each stage. A filter forcing 
1 1 2/1 1 — \un\\ < 10~ 8 was applied to the case A = 2. This cut the path- 
tracking time to less than half a second for N = 17 mesh points. This drastic 
reduction in time for path-tracking as well as data management allowed for 
the confirmation of the existence of two real solutions for up to 100 mesh 
points. Despite the size of the polynomial system when N = 100, each path 
took less than 4 seconds to track from t = 1 to t = 0. A graph of the two 
real solutions for N = 20 mesh points is given in Figure 1. 

3.3 A problem with infinitely many solutions 

It was shown in [3] that the two-point boundary value problem 

y" = -Ay 3 , y(0) = y(l) = 0, (5) 

with A > has infinitely many oscillating real solutions on the interval [0,1]. 
Moreover, the solutions occur in pairs in the sense that —y is a solution 
whenever y is a solution. Hence, together with the trivial solution y = 0, 
we expect always to have an odd number of solutions. That was confirmed 
computationally, as shown in Table 2. Only the case of A = 1 is displayed 
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as all other cases are identical modulo scaling. It may be observed that 
the number of real solutions found by Bertini grows without bound for this 
problem, as the number of mesh points increases. In fact, beyond some 
small value of N, the number of real solutions approximately doubles for 
each subsequent value of N. 



N 


SOLS(iV) 


REAL (TV) 


1 


3 


3 


2 


3 


3 


3 


9 


3 


4 


27 


7 


5 


81 


11 


6 


243 


23 


7 


729 


47 


8 


2187 


91 



Table 2: Solutions of © 



3.4 The Duffing problem 

One representation (see (S]) of the Duffing problem is the two-point bound- 
ary value problem 

y" = -Xsm(y) (6) 

on the interval [0, 1] with y(0) = 0, y(l) = 0, and A > 0. Since our atten- 
tion is restricted to polynomial nonlinearity only, we approximate sin(y) by 
truncating its power series expansion, yielding the problem 

// x / .. y 3 



using two terms or 



using three terms. 

It is known that there are 2k + 1 real solutions to the exact Duffing 
problem © when kir < A < (k + 1) ir. For a given value of A, the 2k + 1 
real solutions include the trivial solution y = and k pairs of solutions 
(yi(x), 2/2(2)) such that y\{x) = —y2(x). Each pair oscillates with a different 
period. As two- and three-term Taylor series truncations for sin(y) do not 
approximate sin(y) well outside of a small neighborhood, the solutions to 
ijTj) and Q may behave quite differently than those of ©• 

Table 3 indicates the number of real solutions found for problems (j7| 
and JHJ) for A = 0.57T, 1.57T, and 2.5n. All solutions have either odd or even 
symmetry about x = i, so we again used the filter ||yi| — |yjv|| < 10~ 8 . The 
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filter was first applied when iV = 4, so the number of real solutions reported 
in each case of Table 3 is the number of real solutions found for N > 5. For 
A = 0.57T and A = 1.57T, there were more real solutions found for (JBJ) than 
predicted for the exact problem ©. However, the computed solutions in 
each case included one pair of solutions that oscillated wildly. These poorly- 
behaved solutions are readily identified by the y'" filter discussed in Section 
2: for N = 25 mesh points, they had residuals four orders of magnitude 
larger than those of the well-behaved solutions. 



A 


f(y) = y- ir 


f(y) = y-ir + m 


f(y) = sin(y) 


0.57T 


1 


3 


1 


1.57T 


1 


5 


3 


2.57T 


1 


5 


5 



Table 3: Number of real solutions for approximations of the Duffing problem. 



3.5 The Bratu problem 

The Bratu problem on the interval [0, 1] has the form 

y" = -\ e y, y(0) = y(l) = 0, (9) 

with A > 0. As in the case of the Duffing problem, we make the right-hand 
side polynomial by truncating the power series expansion of e y , yielding 

y" = -\(l + y+ V ^) (10) 

As discussed in 0, there are two real solutions if A is near zero and no 
real solutions if A is large. The real solutions for small A are symmetric and 
nonnegative. The expected number and properties of the real solutions in 
the cases of A = 0.5 and A = 10 were confirmed, and, as anticipated, 2^ 
total solutions were found in each case for iV = 1, . . . , 15. 

4 Discussion 

A new algorithm for finding the real solutions of a two-point boundary value 
problem has been presented, and several examples have been documented 
under the assumption of polynomial nonlinearity. Furthermore, the use of 
filtering rules to drastically reduce the computational work has been consid- 
ered. In each example presented, the number of real solutions predicted by 
theory has been confirmed computationally, although it was seen that the 
use of filters may effect the number of real solutions discovered. 

There are several variations to the algorithm that could be considered 
in the future. A more detailed analysis of the benefits and drawbacks of the 
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use of filters could be made. Also, it is possible to add extra mesh points at 
the left-hand end or middle of the interval rather than the right. Similarly, 
one new mesh point could be added to each end simultaneously, yielding 
(deg p(y)) 2 starting solutions for each solution from the previous stage. For 
that matter, non-uniform grids could be analyzed with only mild changes to 
the formulation. A similar algorithm could also be developed for systems of 
differential equations. 
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